Nature Computational Science

1Clinical metagenomic sequencing is increasingly used to detect pathogens directly from clinical samples, but sparse and uneven genome recovery caused by low microbial biomass and stochastic fragment loss makes species-level validation difficult. Under these conditions, the same observed coverage may arise either from genuine genome-wide representation or from reads confined to conserved loci shared among related taxa. Current heuristic thresholds based on read counts, breadth or related alignment summaries target the realised coverage rather than the generating process. We show that, under hidden fragment loss, coverage is the wrong inferential object for uniformly valid species-level inference. What remains informative is the spatial dispersion of surviving reads across genomic coordinates. This change of inferential object allows us to ask how much of the candidate genome must be represented in any source capable of reproducing the observed evidence. This yields the Metagenomic Alignment Validation Index (MAVI), a conservative lower bound on the fraction of the candidate genome required to explain the observed read dispersion. MAVI therefore provides an error-controlled, interpretable measure of genomic confidence: the minimum genomic content a competing source would need to share with the candidate to mimic the observed evidence. Evaluated across 643 long-read metagenomic samples, including 280 clinical samples, MAVI reproduced independently supported organism calls, resolved sparse low-coverage ambiguities and exposed workflow-level collapse of stochastic dispersion, without relying on heuristic coverage thresholds.

Inferring the history of effective population size (Ne) from genomic data is fundamental to population genetics, but model-flexible methods have not kept pace with biobank-scale samples. We present Stairway Plot v3, which infers a nonparametric Ne (t) trajectory with bootstrap confidence intervals directly from the site frequency spectrum (SFS). By recasting the optimization as a convex problem and combining inferences across multiple sub-sampled projections, the method scales to 100,000 haploid sequences on commodity CPU hardware without a GPU. Larger samples sharpen the resolution of recent demography, and we apply the method to five human populations including 54,302 Japanese individuals, to our knowledge the largest sample used for nonparametric SFS-based Ne (t) inference.

Single-cell transcriptomic atlases resolve thousands of putative brain cell types, yet common analytical workflows often rely on low-dimensional embeddings that can distort neighborhood relationships and obscure continuous variation between closely related populations. We introduce single-cell Discrete Morse Graph Construction (scDMGC), a deterministic framework based on topological data analysis and discrete Morse theory that extracts a compact graph from neighborhood structure in high-dimensional gene-expression space. Because graph edges are derived from the original k-nearest-neighbor complex, scDMGC provides an interpretable multiscale representation of transcriptomic organization without requiring a low-dimensional embedding. Local maxima represent coherent transcriptional states, saddles quantify their separation, and gradient paths capture intermediate states and continuous transitions. Applied to cortical and hippocampal datasets, scDMGC recovers established inhibitory-neuron classes and elucidates transcriptional gradients. In mouse whole-cortex data, Morse graph structure evaluates cell-type hierarchy and type robustness. Persistence across scales provides a quantitative measure of cell-type identity and discrete versus continuous relationships. Finally, in Alzheimers disease single-nucleus data, scDMGC distinguishes changes in cell-type abundance from disease-associated shifts in transcriptional state. scDMGC therefore provides a multiscale framework for quantifying whether transcriptional populations form persistent cell types, how strongly external annotations are supported by intrinsic data structure, how cell types relate to one another, and where discrete organization transitions into continuous variation.

Mendelian randomization (MR) uses genetic variants as instrumental variables for causal inference, but its validity is often undermined by horizontal pleiotropy and linkage disequilibrium (LD) among instruments, which current methods typically address by discarding correlated or outlying variants. Here we present ECC-MR, an MR framework that brings ideas from error-correcting codes into causal inference: a parity-check matrix constructed from LD structure places the causal signal in its null space, so that pleiotropic effects become a decodable sparse error syndrome rather than a nuisance. We establish theoretical guarantees for consistency and stable recovery under sparse pleiotropy, and show in simulations that ECC-MR outperforms existing methods in bias, efficiency and Type I error control, while flagging regimes where its assumptions fail. Applied to three exposure-outcome pairs spanning weak to strong pleiotropic burden, ECC-MR yields concordant causal estimates while retaining all correlated instruments, and its decoded pleiotropic variants recover known hotspots and drug targets, corroborated by colocalisation and phenome-wide association analyses. ECC-MR thus turns LD and pleiotropy from liabilities into resources for robust causal inference and mechanistic discovery.

Focal cortical dysplasia (FCD) is among the leading causes of drug-resistant focal epilepsy and can be challenging to detect on routine MRI due to the subtle and heterogeneous appearance of lesions. Motivated by the scarcity of clinically representative cohorts, we investigate the feasibility of simulating FCD-like lesions in healthy MRI using a composition of hand-crafted image transformations. We develop SynthFCD, a phenomenological FCD simulator, and create 2,164 synthetic examples from healthy controls in the heterogeneous FOMO300k dataset. We then pretrain deep learning (DL) segmentation models on synthetic lesions and fine-tune them on real FCD cases from the University Hospital Bonn (UHB) cohort. Comprehensive evaluations show consistent improvements across simulator presets, with detection-rate gains of up to 18.5 percentage points over training from scratch and 7.4 points over self-supervised pretraining, while achieving state-of-the-art segmentation performance (0.366 DSC) on the challenging held-out UHB benchmark. These findings suggest that simulation-based synthetic data generation can bridge existing radiological knowledge and advances in DL to provide meaningful supervisory signals that could ultimately aid the detection of underrepresented and diagnostically challenging cases.

We present a lossless compression algorithm for indexing a protein database while supporting fast taxonomic classification in the method Centrifuger. The algorithm is a new scheme of the previously proposed run-block compression algorithm to reduce the size of the Ferragina-Manzini (FM) index, and it scales better with alphabet size than the original. On the RefSeq prokaryotic and viral protein sequences, Centrifuger reduces the memory footprint by over a third compared to the method Kaiju that builds on a plain FM-index, while having comparable running time. Furthermore, the compressed FM-index is lossless and can locate matches of arbitrary length, which helps Centrifuger achieve greater accuracy than Kraken2, a k-mer-based taxonomic classification method. We leverage the computational efficiency of Centrifuger to create an index of size 182 GB for classifying the reads against the full nr database that contains about 250 billion amino acid characters. Using this index, Centrifuger reveals different SARS-CoV-2 infection states and viral transcriptome profiles across human cell types from single-cell RNA-seq data.

Generating diverse, all-atom protein conformations remains challenging, especially for recovering alternative functional states and modeling conformational ensembles. We present trFlow, an ultrafast approach for generating all-atom conformations using geometry-guided flow matching. trFlow leverages three core innovations: (1) a hybrid single- and multi-step sampling scheme that balances accuracy and speed; (2) a decoupled representation-generation architecture that maximizes computational efficiency; and (3) fine-tuning on NMR ensembles to better capture conformational heterogeneity. Extensive benchmarking across dual-conformation proteins, NMR ensembles, and a CASP16 dynamic target shows that trFlow outperforms competing methods in recovering challenging states and reproducing experimental heterogeneity, while being an order of magnitude faster than existing all-atom approaches. These results establish trFlow as a practical and scalable framework for large-scale generation of all-atom conformational ensembles.

Transcranial direct current stimulation (tDCS) is a promising non-invasive therapy for drug-resistant epilepsy, but the montage optimization pipelines used to personalize it often rely on biophysical head models that capture only the induced electric field and ignore the network dynamics through which stimulation influences epileptogenic activity. We aimed to develop and evaluate a montage optimization framework that incorporates individualized network dynamics. We introduce a framework based on a digital twin, or neurotwin, a personalized replica of a patient's brain that couples a biophysical head model of the stimulation-induced field with a personalized whole-brain model (WBM) of seizure dynamics constrained by structural connectivity and intracranial recordings. We applied it to produce optimal montages in 12 patients with drug-resistant epilepsy, compared them with conventional biophysics-based optimization, and assessed the framework retrospectively against the clinical outcomes of six treated patients. Because the neurotwin optimization suppresses simulated spread by construction, the informative result is the divergence between pipelines: the network-informed and field-based optimizations produced substantially different montages despite delivering comparable inhibitory fields at the epileptogenic zone, and montages with nearly identical field distributions could yield markedly different network outcomes within the model. A minimum-replacement-set analysis traced these differences to field changes at a few propagation-zone and off-target parcels rather than to the overall field. The parcels whose inhibition most reduced spread tended to occupy highly connected (hub) positions in the structural connectome, a robust but partial association indicating that network structure shapes, but does not fully determine, where stimulation is most effective. In the retrospective analysis, model-predicted spread reduction tracked clinical seizure-frequency change in the expected direction, though not significantly (n=6). We hypothesize that accounting for individualized network dynamics may improve seizure control, a hypothesis now being tested prospectively (NCT06334952).

Large-scale biophysically-detailed models of the brain are instrumental in unraveling how cellular, synaptic, and circuit-level mechanisms contribute to neural dynamics. Such models are built around large populations of single-neuron models whose morphologies must preserve physiological variability to reproduce integrative properties and firing patterns observed experimentally. Experimental reconstruction of morphologies is a time-consuming and resource-intensive process, unlikely to supply them at the required scale. Several data-driven methods for morphology synthesis have been developed, which nonetheless rely on extensive datasets for parameter estimation. To address this limitation, we developed NeuWalk, to our knowledge the first public Python library for morphology synthesis based on branching-and-annihilating, biased-and-correlated random walks. It estimates branching and annihilation rates from Sholl intersections and, optionally, bifurcation counts, while also ensuring their replication in synthetic morphologies; user-specified combinations of biases, implementing somatofugality, self-avoidance, spatial competition, and branching-angle control, sculpt morphological patterns characteristic of each neuron class. We then validated the method on mitral and middle tufted cells of the olfactory bulb, semilunar and pyramidal cells of the anterior piriform cortex, and neocortical pyramidal neurons: means and variances of Sholl intersections, bifurcation counts, and total dendritic length were statistically indistinguishable between synthetic and natural morphologies, with the single exception of its variance in apical dendrites of piriform cortex pyramidal neurons. As Sholl intersections and bifurcation counts are commonly reported in the literature, and biases can be inferred from published images or qualitative descriptions, NeuWalk requires no immediate access to experimental reconstructions, striking a trade-off between data-driven parameterization and empiricism without sacrificing morphological accuracy.

Summary: Ensembl VEP is the de facto reference for variant consequence annotation, but its Perl implementation limits throughput, and the open reimplementation benchmarked here is validated on a set too small to bound its disagreement with it. We present vep-rs, a Rust reimplementation of Ensembl VEP release 115.2, audited over 260,638,766 VEP consequence tuples, each a (location, allele, feature) key carrying a Sequence Ontology term set, from six variant datasets (ClinVar, gnomAD chr21, and 1000 Genomes chr21 across GRCh37 and GRCh38). On SNPs and indels vep-rs emits exactly VEP's tuple count on all six and attains F1 [≥] 0.999974, differing from VEP on 3,169 tuples, at most one in 38,000 per dataset; with two documented VEP defect shapes set aside, F1 rounds to 1.000000 on all six. Per consequence class, 29 of the 30 Sequence Ontology terms carrying at least 1,000 tuples exceed F1 0.99, and the lowest, 0.540 on start_retained_variant, is entirely a VEP defect. On structural variants, F1 is 0.975395 (GRCh37) and 0.909998 (GRCh38). Timed on 20 independent machines per engine and architecture, vep-rs is faster than VEP by 101x-284x on ARM (geometric mean 176x) and 78.3x-192x on x86 (geometric mean 135x). Five divergence classes are documented, four of them defects in VEP release 115 that vep-rs does not reproduce: a consequence call that contradicts itself, a splice-region term dropped from a variant that lies in a splice region, annotation against a chromosome the variant does not lie on, and output that depends on input batch composition. A pipeline keyed on VEP consequence terms can therefore adopt vep-rs against a measured, per-term bound, shed VEP's contradictory and batch-dependent calls, and annotate population-scale inputs at 78.3 to 284 times VEP's throughput, and per-file structural-variant sets at 22.5 to 52.9 times. Availability and implementation: vep-rs, its harness, both comparators, and every per-machine wall time are released under Apache-2.0 at https://github.com/natera-open-source/vep-rs (release v0.1.0, archived at https://doi.org/10.5281/zenodo.22837897). Supplementary information: Supplementary data are available at Bioinformatics online.

Flow cytometry is central to clinical hematology, generating extensive archives comprising billions of single-cell profiles linked to expert-derived diagnostic summaries. While these records represent a largely untapped repository of collective expert knowledge, the underlying single-cell annotations required for supervised machine learning are typically not recorded in a standardized, accessible format. We present IMPGAT (implicit gating), a weakly supervised multiple-instance learning framework that reconstructs single-cell gating strategies directly from existing sample-level cell-type frequency distributions. The model utilizes an induced set attention mechanism to map single-cell marker profiles to specific cell populations, requiring neither explicit cell-level labels nor unsupervised clustering during training. We evaluate the framework on acute lymphoblastic leukemia (ALL) and acute myeloid leukemia (AML) cohorts, totaling 711 measurements from 268 individuals and 240 million single cells. In ALL, the model's quantification of cell-type abundances shows no statistically significant difference from human inter-expert variability. In AML, the framework achieves blast quantification on par with a recently established supervised state-of-the-art model that relies on explicit single-cell annotations for training. Error analysis shows that mismatches between model predictions and expert gating primarily involve expert-labeled non-cellular events, reflecting the non-incorporation of any explicit physical prior knowledge (such as the geometric relationships traditionally used to identify doublets and debris in scatter cytograms) in model training. Conversely, model certainty correlates with phenotypic borderline cases where manual bisection is inherently difficult. Our approach enables the systematic utilization of retrospective clinical records for the development of automated, standardized gating logic and is applicable to broader high-throughput contexts, such as large-scale immune monitoring or phenotypic drug screening.

Mapping the multi-mutant fitness landscape is vital to protein engineering, but is challenging due to the vast combinatorial sequence space awaiting exploration. A central difficulty lies in the accurate and efficient modeling of non-additive epistatic effects among individual mutations, which partially arise from the physical inter-residue interactions prescribed by the protein structure. Existing fitness predictors usually perform well on single mutants but become less powerful for higher-order mutants, due to the lack of explicitly considering the relationship between sequence, structure and function of the target protein. Here, we present an end-to-end framework named Cerebra-Epistasis, which couples a single-sequence structure predictor that explicitly endows the structure awareness beyond the conventional sequence-fitness mapping with a downstream fitness prediction network that deliberately models the non-linear epistatic effects beyond the traditional additive terms. When evaluated across diverse assays, Cerebra-Epistasis outperforms the other state-of-the-art baselines in multi-mutant fitness prediction, with enhanced advantage over increasing mutation orders. Moreover, our special design on the epistasis modeling allows reliable extrapolation from low-order mutant data to unseen higher-order combinations, enabling one-shot inference of the overall mutation atlas from the starting sequence, a benefit that supposedly introduces three orders of magnitude acceleration in the landscape-scale prediction.

Motif scaffolding builds a protein around a prescribed structural motif -- an active site, a binding epitope, a metal site -- so that the motif is presented in its required geometry. The task is usually posed with the motif pinned down: its position along the sequence is given, and several rigid fragments inherit their mutual pose from a native structure. Many applications supply neither. We present NomadFlow, unindexed multi-motif scaffolding on all-atom SE(3) flow matching: the motif enters as a condition, a group projection preserves the internal geometry of each group exactly, and the groups move freely relative to one another until a set time, after which they act only as a condition. We assemble a preliminary benchmark for this setting, varying the granularity of the partition and the structural context around the motif, and obtain unique successes on 19 of 21 tasks. The trained model proposes a diversity of valid mutual poses, and on several tasks the native arrangement is among them. The shape of that solution space, however, varies strongly with the geometry of the condition: across tasks the successful layouts form a wide cloud, collapse onto the native pose, or contract to a connected low-dimensional family -- a continuum of poses rather than a set of isolated solutions. Multi-motif methods are therefore better compared by the space of layouts they propose and the successful part of it they retain than by yield or fold diversity alone. Code and weights: https://github.com/ipermyakov/nomadflow.

Local ancestry inference (LAI) reveals population compositions of an individuals chromo-somes, which empowers population genetics and genetic association studies of diverse populations. The recent availability of large genomic datasets from biobanks offers unprecedented opportunities for LAI. However, such datasets are typically only accessible from cloud-based platforms, and the lack of efficient and cost-effective tooling makes LAI prohibitive for resource-constrained researchers. Here, we develop Recomb-Mix Cloud, a cloud-ready pipeline for LAI based on Recomb-Mix, a recently introduced efficient LAI method. Recomb-Mix focuses on individual variants while maintaining a lower computational burden than other LAI methods focusing on genomic windows. Recomb-Mix Cloud offers a cloud-friendly complete pipeline from panel sample selection, marker alignment, compact reference panel construction, Recomb-Mix threading, and generating LAI calls in common data formats. Recomb-Mix Cloud also offers a set of pre-computed compact population reference panels constructed based on data from the Human Genome Diversity Project and the 1000 Genomes Project. Our results showed Recomb-Mix is about one and two orders of magnitude cheaper and faster than FLARE and RFMix, two popular LAI methods, while offering largely consistent LAI calls on the All of Us Researcher Workbench. The estimated overall compute cost per All of Us sample for whole-genome LAI analysis using Recomb-Mix Cloud with a pre-computed reference panel is just $0.01, compared with $0.43 for FLARE and $4.58 for RFMix. Therefore, Recomb-Mix Cloud is well-positioned to democratize local ancestry-based genetics research on modern biobank-sized data. The Recomb-Mix Cloud pipeline and the compact reference panels are available at https://github.com/ucfcbb/Recomb-Mix.

High-accuracy, high-resolution benchmark data are fundamental to improving the accuracy and generalizability of sequencing technologies and associated algorithms. Most existing benchmarks for DNA sequencing rely on mapping to a single linear reference, omitting inaccessible and highly repetitive complex genomic regions. Here, we present the Genome in a Bottle (GIAB) Consortium's HG002 v5.0q benchmark set, a more comprehensive variant benchmark derived from a fully phased, highly curated, and accurate diploid assembly of HG002. Our approach leverages recent advances in long-read sequencing and assembly to provide confident small- and structural-variant calls across previously inaccessible regions, including segmental duplications, homopolymers, tandem repeats, and complex immune loci. The v5.0 benchmark significantly expands the callable genome relative to v4.2.1, adding x200 Mbp of sequence and increasing variant detection by 18 % for small variants and 3-fold for structural variants compared to previous benchmarks. The benchmark regions include more variants in repetitive regions, including 1.2 x segmental duplications, 1.8 x tandem repeats, and 2.0 x homopolymers, compared to v4.2.1. Separate and combined small- and structural-variant benchmarks are available for three versions of the human reference genome (GRCh37, GRCh38, and CHM13). We validate the reliability of this benchmark set by comparing multiple sequencing technologies and variant calling pipelines and curating differences. The v5.0q benchmark provides a critical resource for the development and validation of next-generation sequencing technologies and analysis methods.

Estimating population-level dependence between brain regions is challenging when the number of regions greatly exceeds the number of individuals. We introduce an information-borrowing generalized Bayesian method for high-dimensional covariance estimation from region-level brain measurements. A divide-and-conquer parameterization separates marginal variances from correlation calibrations: an empirical \emph{pilot} correlation matrix--by default, obtained by standardizing a sample covariance matrix (SCM)--provides a reference structure that is adaptively calibrated through region-specific parameters shared across pairwise correlations, reducing the number of free parameters from $\mathcal{O}(d^2)$ to $\mathcal{O}(d)$ ($d$ is the number of regions) while guaranteeing a positive-definite covariance through the inclusion of a ridge term. Horseshoe priors on the calibration parameters enable selective borrowing and shrinkage of pilot information. Under suitable conditions on pilot accuracy, we establish theoretical guarantees for covariance estimation while retaining scalability in ultra-high-dimensional settings via a low-rank likelihood representation and a Woodbury-accelerated Hamiltonian Monte Carlo algorithm. Simulations showed lower estimation error than the SCM, bagging and iteratively refined principal orthogonal complement thresholding (POETRY). Applications to a semi-synthetic cortical-thickness and ADHD resting-state fMRI data revealed significant group differences in the default-mode--visual network and somatomotor network--limbic network blocks after Benjamini--Hochberg correction.

Generative models of protein structures have accelerated scientific discovery. However, the leading structure-prediction models are trained primarily on crystal structures, limiting their ability to sample the conformations of intrinsically disordered proteins (IDPs). All-atom molecular dynamics (MD) simulations can probe the diversity of IDP conformational ensembles, but these simulations are slow and costly. Leveraging the speed of generative models and the accuracy of MD, we develop Zephyr, an SE(3) equivariant diffusion model trained on simulated IDP ensembles. The model generates the structure of a single protein sequence in two phases: a coarse-grained structure is first sampled and heavy atoms are subsequently denoised onto this structure. We evaluate our method on two fronts: ensemble accuracy and sampling efficiency. Across five test systems, the model reproduces key geometric and energetic metrics of the MD ensembles, while sampling at least 63 times faster.

Cellular diversity is required for the proper function of many tissues, including the retina, with its ~130 cell types.1-3 Previous lineage studies suggested that retinal progenitor cells (RPCs) generate cell diversity through biased cell divisions.4-9 However, the extent to which each progenitor cell contributes to diversity, and the overarching rules governing the lineages of the collection of progenitor cells, are unknown. While classic lineage tracing approaches including retroviral labeling4,6,10, dye injection11,12, and time-lapse microscopy13,14 were used to make the foundational insights regarding retinal lineages, they cannot be used to interrogate multigenerational lineages at scale. Here we describe the creation of an evolvable lineage tracing tool, "SCRIBBLE" for Sequential Combinatorial Recorder for Iterative Barcode-Based Lineage Evolution. SCRIBBLE uses prime editor to iteratively edit a barcode sequence while preserving the order of genesis of clonally related progeny. Through a series of rationally designed high-throughput screens and design optimizations, we identified novel sequences that support a high level of barcode editing and information content. Evaluation of SCRIBBLE in embryonic mouse retina recapitulated known lineages in retinal development, demonstrating that SCRIBBLE is capable of high-resolution lineage tracing in a complex tissue.

Online identification of nonlinear dynamical systems faces two challenges. Noisy, partial observations make the latent state difficult to estimate. Local attractors can also keep passive trajectories in regions that reveal little about the unknown dynamics. Both challenges commonly arise in closed-loop neural experiments where stimulation must be chosen from noisy spike counts, while neural activity may remain near a stable state unless stimulation drives it elsewhere. We address this problem by actively selecting inputs to collect informative observations during the experiment. We introduce PALDI (Planning-based Active Latent Dynamics Identification), which combines extended-Kalman state prediction, single-step Laplace state updates, blockwise Gaussian parameter updates, and finite-horizon planning. Its parameter expected-information-gain surrogate combines dynamics sensitivity with state Fisher information and discounts apparent parameter information when latent states are uncertain. We compare PALDI with nonadaptive inputs, a myopic policy, and two active baselines: the one-step D-optimal method FLEX and receding-horizon curiosity with uncertainty sampling. On a bistable Wilson-Cowan model of excitatory and inhibitory populations, PALDI approaches the true-model reference and remains strongest as the observation signal-to-noise ratio decreases or the readout geometry becomes biased. We then test identification in a spiking decision network where we use a reduced dynamics and observation model to approximate the data-generating process. PALDI reaches a median held-out prediction R2 of 0.90, close to a fitted reference, and its model gives the highest change-of-mind success among identification policies. Information-guided finite-horizon planning is a practical strategy for identifying neural circuit dynamics with closed-loop stimulation.

Whole brain tractography contains rich trajectory-level information about white matter organization, yet conventional structural connectomes capture only part of it. We introduce Structural Network Coupling (SNC), a network representation enabled by SILICA (Streamline Independent Component Analysis). SILICA decomposes tractography into structural source components and their streamline-specific loading profiles, while SNC characterizes relationships among these components through their shared trajectory expression. The resulting component-level network retains direct correspondence to the underlying white-matter pathways and supports network visualization and graph-based analysis. SNC produced anatomically interpretable coupling among trajectory-resolved white-matter components and revealed a highly organized group architecture that independently recapitulated commissural, projection, and association systems. It also showed strong hemispheric segregation and structured cross-system coupling, and was highly stable across split-half samples. SNC extends SILICA from trajectory decomposition to network-level organization and provides a flexible framework for future structural and multimodal network analysis.